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ABSTRACT 


Two new approaches are proposed and developed for making time and space- 
dependent, quantitative short-term forecasts of lightning threat, and a blend of these 
approaches is devised that capitalizes on the strengths of each. The new methods are dis- 
tinctive in that they are based entirely on the ice-phase hydrometeor fields generated by 
regional cloud-resolving numerical simulations, such as those produced by the WRF model. 
These methods are justified by established observational evidence linking aspects of the 
precipitating ice hydrometeor fields to total flash rates. The methods are straightforward 
and easy to implement, and offer an effective near-term alternative to the incorporation of 
complex and costly cloud electrification schemes into numerical models. 

One method is based on upward fluxes of precipitating ice hydrometeors in the mixed 
phase region at the -15° C level, while the second method is based on the vertically in- 
tegrated amounts of ice hydrometeors in each model grid column. Each method can be 
calibrated by comparing domainwide statistics of the peak values of simulated flash rate 
proxy fields against domainwide peak total lightning flash rate density data from observa- 
tions. Tests show that the first method is able to capture much of the temporal variability 
of the lightning threat, while the second method does a better job of depicting the areal 
coverage of the threat. Our blended solution is designed to retain most of the temporal 
sensitivity of the first method, while adding the improved spatial coverage of the second. 

Simulations of selected diverse North Alabama cases show that WRF can distinguish 
the general character of most convective events, and that our methods show promise as a 
means of generating quantitatively realistic fields of lightning threat. However, because the 
models tend to have more difficulty in predicting the instantaneous placement of storms, 
forecasts of the detailed location of the lightning threat based on single simulations can 
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be in error. Although these model shortcomings presently limit the precision of lightning 
threat forecasts from individual runs of current generation models, the techniques proposed 
herein should continue to be applicable as newer and more accurate physically-based model 
versions, physical parameterizations, initialization techniques and ensembles of forecasts 
become available. 
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1. Introduction 


The threat from lightning in convective storms is a significant source of concern for 
public safety and a wide range of weather sensitive operations. In the United States 
alone, lightning is responsible for nearly 1000 deaths and injuries each year, with damages 
exceeding SI billion (Curran' et al. 2000). As a consequence, improved short-term (0-12h) 
forecasts of lightning threat are of interest to the National Weather Service (NWS) and 
other forecasting organizations (Darden et al. 2006). Although cloud-to-ground (CG) 
lightning is of obvious importance to ground-based operations, the total lightning activity 
(intracloud and CG) has been shown to be useful as an indicator of impending severe or 
high-impact weather (Goodman et al. 2000; Williams et al. 1999; Gatlin 2007). With the 
planned launch in 2014 of a Geostationary Lightning Mapper (GLM) aboard the GOES- 
R series of Geostationary Operational Environmental Satellites, the exploitation of new 
applications derived from total lightning measurements will continue to gain interest. 

Previous efforts to forecast lightning threat have been based largely on the climatolog- 
ical connections between thunderstorm occurrence and certain parameters of the pre-storm 
environment. Bright et al. (2004), for example, utilized midlevel layer estimates of convec- 
tive available potential energy (CAPE) to infer likely convective storm upward velocities 
aloft in their construction of lightning probability maps. However, their approach was 
used only for finding regions favorable for thunderstorm formation, and could not be used 
to make quantitative estimates of lightning flash rate within storms. Regression-based 
methods proposed by Mazany et al. (2002) and Burrows et al. (2005) are similarly lim- 
ited. While Bothwell (2005) and Shafer and Fuelberg (2007) also use statistical methods 
to relate lightning occurrence to aspects of the observed storm environment, they further 
attempt to classify the amounts of lightning in subjectively defined categories. All these 
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previous research efforts have focused on the threat from cloud-to-ground lightning. 

In recent years, global-scale observing systems such as the Tropical Rainfall Measuring 
Mission (TRMM) satellite, with its Lightning Imaging Sensor (LIS) and Precipitation 
Radar (PR), have allowed investigators to document robust relationships between lightning 
flash rates and the hydrometeor microphysical properties of storms. Cecil et al. (2005) 
showed that the flash rates of storms were closely tied to the amount of radar reflectivity 
aloft in the ice-phase regions of the storms. Meanwhile, Petersen et al. (2005) found 
a strong relationship between the vertically integrated large precipitating ice content in 
storms and their lightning flash rates. Fluxes of ice species such as graupel in the mixed 
phase region near -15°C have also been shown to be related to storm flash rates (Petersen et 
al. 2005; Deierling and Petersen, 2008). These regime-independent globally-based results 
are consistent with the more limited regional observations reported by Shackford (1960), 
Goodman et al. (1988; 2005), and others. 

Most cloud-resolving numerical models now have the capability of computing fields of 
mixing ratios of multiple species of hydrometeors, including several important ice-phase 
species. Thus, cloud-resolving models, such as the Weather Research and Forecasting 
(WRF; Skamarock et al. 2005) model, can now provide time- and space-dependent sim- 
ulated fields of parameters known to be associated with lightning flash rate. Indeed, 
cloud-resolving forecast simulations with these capabilities are now becoming common- 
place. However, in recognition of the rapid growth of errors in models at convective scales, 
some recent simulation investigations have even begun to explore the use of cloud-resolving 
model ensembles (see, e.g., Kong et al. 2006). 

In this paper, we develop methodologies to demonstrate how regional cloud-resolving 
forecast simulations can be exploited to create quantitatively calibrated, time-dependent 
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and specific short-term forecast maps of evolving lightning flash rate density fields in 
convective environments. Given high-resolution output from a suitable numerical model, 
our prototype methods yield lightning forecast products that are straightforward, while 
avoiding the added expense and complexity of incorporating explicit cloud electrification 
algorithms into the models (see, e.g., Helsdon et al 1987, Helsdon et al 1992, MacGorman 
et al 2001, Mansell et al. 2002, Kuhlman et al 2006). For this prototype research, 
we restrict our efforts to analysis of a selection of deterministic simulations only. The 
application of our methods to convective model ensembles is considered beyond the scope 
of this work, and is relegated to future research. 

Two distinct approaches to forecasting the lightning flash rate density field are devised 
and discussed herein, along with a blended version that attempts to capitalize on the 
strengths of each. One uses the prognosed field of upward vertical velocity multiplied by 
graupel mixing ratio at the -15°C level, hereinafter called “graupel flux,” even though mass 
units are absent because air density is not used in its evaluation. The second approach 
exploits the relationship between storm total flash rate and total volumetric amount of 
precipitating ice. For use on a. model grid, we may express this latter relationship in 
terms of the vertically integrated ice content in each grid column. Other approaches 
based on regression analysis of satellite-based lightning flash rates against radar reflectivity 
profiles in storms (see Cecil et al 2005), are also under consideration, but require further 
analyses of the TRMM PR-LIS databases and additional simulation cases for algorithm 
development. Both of the techniques reported herein can be empirically calibrated for any 
given cloud-resolving model configuration, using successful simulations of a few diverse 
cases for which ground-truth total lightning flash rate densities are known. Because our 
blended approach is simply a weighted average of methods 1 and 2, it maintains the correct 
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calibration for peak flash rate densities, while retaining significant temporal sensitivity and 
ensuring improved capabilities in matching threat areal coverage to the observations. 

Results of WRF simulations of various types of Alabama storm systems demonstrate 
the model’s ability to discriminate between convective events of differing depth and inten- 
sity, although with the 2 km grid mesh employed here, the model often produces storms 
that are somewhat weaker and exhibit more areal coverage than those observed. The 
model’s weak bias is compensated for by our flash rate calibration methods, while the 
model’s tendency to produce too much areal coverage of convection does not appear to 
extend to the deeper convection that is associated with lightning. However, the model also 
tends to mislocate storms spatially and temporally, and it is not easy to compensate for 
such phase errors. Because of the limitations of the WRF configuration used here, and of 
the small sample of cases studied, our calibrations and results here should be regarded as 
exploratory. Nevertheless, the flash rate prognostic techniques described here show much 
promise, and can readily be applied to improved future versions of cloud models. 

Section 2 explains the data and analysis techniques employed in this study, and 
presents the basis for our calibrated lightning threat methods. Section 3 applies cali- 
bration techniques to our full sample of simulations, thereby obtaining detailed formulae 
for lightning threat computations, and applies these formulae to two contrasting WRF 
model case studies. In section 4, we present our blended solution for WRF-based lightning 
threat, and show how it works for our two selected case studies. In section 5, we give a 
discussion of our results, and in section 6, we summarize our findings and discuss future 
research needs. 

2. Methodology 
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a. Numerical simulations 


Our numerical simulations are conducted using the WRF model, version 2.1.2, on 
a native 2-km horizontal mesh having 51 vertical sigma levels. For this model mesh, 
convection is simulated explicitly; no cumulus parameterization is needed. Other model 
physical parameterization choices are detailed in Table 1. Of particular importance is 
the use of the WRF Single-Moment Six Species (WSM-6) microphysics package. In this 
simplified scheme, water substance can assume only the six forms of vapor, cloud water, 
rain, cloud ice, snow and either graupel or hail. In our model runs, we have opted to use 
graupel, with a density of 300 kg m -3 . 

All model runs are initialized using AWIP212 National Centers for Environmental 
Prediction (NCEP) Eta Data Assimilation System (EDAS) analyses, supplemented with 
surface and aircraft observations, and also data from the national network of WSR-88D 
Doppler weather radars. The radial velocity fields from the WSR-88Ds were always used, 
and reflectivity fields were also tested in some cases where precipitating convection was 
already occurring at the model start time. Three-hourly forecast fields from NCEP’s Eta 
model were used to obtain updates of model lateral boundary conditions. 

WRF integrations were conducted for periods lasting from 6 to 12 h, with 25 time 
periods saved for analysis. The save interval was 15 min for the 6-h forecast runs. For the 
12-h runs, only the last 6-h of the forecast output was saved. This was done because, for the 
12-h runs, most of the significant convection occurred more than 6-h into the simulation. 
One simulation, for the 10 December 2004 cold season hailstorm case discussed herein, 
was run for 8 h, with output saved at 20 min intervals. Model data were saved as binary 
files interpolated to a latitude-longitude grid centered near Huntsville, Alabama, and were 
analyzed using the Grid Analysis and Display (GrADS) software. The mesh spacing for 
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the interpolated grid was approximately 0.009 degrees in both directions, which equates 
to 1 km spacing along meridians, approximately twice as fine as the native mesh used in 
the WRF simulations. The “grid boxes” used in our analyses are therefore slightly smaller 
than 1 km 2 . The grid mesh contained 447 elements in the zonal direction, and 335 in 
the meridional direction. Model data were also interpolated in the vertical to yield fields 
with 500 m vertical mesh spacing, with the lowest level located at 250 m above sea level. 
On our native 2 km grid, the high terrain of the Appalachian Mountains intrudes several 
vertical levels into our simulation domain. As a result, all atmospheric data from points 
located beneath the earth’s surface were treated as missing, and were excluded from Our 
calculations of atmospheric quantities. 

Several diverse case studies were selected from the North Alabama region, based on 
the availability of both model initialization and boundary data and also of “ground-truth” 
lightning data from the North Alabama Lightning Mapping Array (LMA; Goodman et al. 
2005), described briefly below. Results from the simulations showed varying degrees of 
success. Simulations that we consider as successful featured deep convection of the proper 
general intensity, with subjectively identified main reflectivity features that appeared to 
resemble those of the observed storms, and which occurred in the right general area (spa- 
tial postioning errors 0(100 km)), at the right time (temporal errors 0(2 h)). In two other 
cases, the WRF model failed either to initiate or sustain convection of the proper intensity, 
and thus did not provide data that could be used to characterize or calibrate the lightning 
threat associated with the convection. For such cases, it was found that even model runs 
initialized with radar reflectivity fields failed to generate sustained intense convection. In- 
sofar as the purpose of this research is to demonstrate methods of using cloud-resolving 
model output as a means of obtaining quantitatively calibrated lightning forecasts - and 
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not to pursue a validation study of any specific model - we consider here only the successful 
WRF simulations. For the purposes of documenting our flash rate density forecast calibra- 
tion methods, we will present results using data from all seven of our successful simulations. 
For the purposes of showing graphically the character of individual forecasts, however, we 
will focus on results from two strongly contrasting cases, one featuring a tornadic supercell 
ahead of a frontal squall line on 30 March 2002, and another featuring smaller and weaker 
cold-season hailstorms on 10 December 2004. The 30 March 2002 event represents a high 
flash-rate case, while the 10 December 2004 case represents the low flash-rate end of the 
spectrum. 

b. Lightning observations 

Our ground-truth lightning dataset consists of flash analyses derived from data col- 
lected by the North Alabama Lightning Mapping Array (LMA; see e.g., Krehbiel et al 
2000). Located in the greater Huntsville, Alabama area, this array consists of 10 sen- 
sors that detect VHF pulses from lightning channel segments. After the sensors record 
the times of received pulses, supplied software (Version 6.2) cross-correlates them and 
identifies which ones likely correspond to individual physical radiation events. Using this 
information, the actual time and three-dimensional location of the physical source events 
comprising lightning flash segments can be estimated (Hamlin, 2004). A flash clustering 
algorithm is then applied to these chronologically-ordered source data to segregate the 
large numbers of individual radiation sources into discrete lightning flashes. 

The flash algorithm employs time and space proximity criteria to identify which 
sources are likely to be part of any given flash. Sources are assumed to be part of the 
same flash if they occur less than 0.3 s apart in time, and also satisfy a spatial separa- 
tion restriction. This spatial separation criterion is spatially inhomogeneous, to reflect the 
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differing range and azimuth dependences of the LMA location error, and the range depen- 
dence of its detection efficiency characteristics. In addition, consideration is also given to 
the maximum likely physical spacings between radiating segments of lightning channels, 
but these spacings are treated as spatially homogeneous, at least in the horizontal. Previ- 
ous research (Boccippio et al. 2000; Thomas et al. 2004; Koshak et al. 2004), has shown 
that the source range and altitude uncertainty increase quadratically with range, while the 
azimuthal. error is approximately independent of range. As with the range location uncer- 
tainty, the LMA detection efficiency decreases rapidly with range. The physical spacing 
between radiating segments of lightning channels, meanwhile, can exceed 1 km. 

Thus, consideration of the above factors, along with study of the radial stretching 
of LMA-derived storm footprints with range, suggests that it is appropriate to assume 
that the maximum allowable total range location difference between successive radiation 
sources that belong to the same flash increases with the square of range from the network 
center, such that the range difference is 10 km at a range of 100 km, and 40 km at a 
range of 200 km. Maximum azimuthal location uncertainties appear to be bounded by 
0.05 rad, which is equivalent to 5 km at a range of 100 km. Thus, pairs of sources that 
satisfy our time difference criterion and also are spaced at ranges that differ by less than 
our empirical range uncertainty and also at azimuths that differ by less than our empirical 
azimuthal uncertainty are considered to be part of the same lightning flash. Note that 
our empirically-derived source-linkage spatial criteria for points in flashes are considerably 
larger than the simple root-mean-square range and azimuthal errors of individual LMA 
points documented in Thomas et al. (2004). This is due to the need to correctly group 
most, if not all, the sources belonging to individual flashes, which requires the use of 
proximity criteria that are much larger than the root-mean-square point location errors of 
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individual sources. The proximity criteria are also made larger by the range dependence 
of detection efficiency falloff, and to the inherent variability of the radiative structure of 
lightning channels. Although our proximity criteria appear generous, tests confirm that 
they are approximately the smallest that can be chosen without causing spurious increases 

in flash counts. Furthermore, comparisons with an interactive algorithm developed by New 

) 

Mexico Tech (see Thomas et al 2003) show that our algorithm’s flash counts are within 
about 5% of theirs, with our algorithm being slightly more conservative. 

Because of the rapid increase in range location uncertainty with range from the net- 
work center, the ability to map fine-scale spatial structure of flashes diminishes with in- 
creasing range from Huntsville. This, however, is not an issue with the present research, 
where only peak flash count densities are used, and the storms producing these peak val- 
ues were within 200 km range of the network center. Because of the nature of our flash 
clustering algorithm, we believe any errors in our flash counts would be manifest as only 
negligibly small undercounts. These small undercounts are at least partly offset by the 
possible splitting of some flashes that straddle the beginning times of our 5-min lightning 
analysis periods. 

One issue that complicates interpretation of the derived flash rates is the occurrence of 
what appear to be single-source flashes. These ’’singleton” events may be due to errors in 
correlating weak raw signals identified as being parts of an individual meteorological source, 
or to errors involving mixing of nonmeteorological pulse detections with meteorological 
ones (Thomas et al. 2004). However, because many singletons are observed to cluster 
near meteorological targets such as storms, they may also correspond to weak electrical 
discharges that are not directly associated with stronger, more easily recognized flash 
events. Because their validity as true independent lightning events is still open to question, 
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and because they stand out from what appears to be a true background continuum in the 
spectrum of source numbers per flash, we exclude all singletons from the analyses used in 
this paper. Thus, we consider valid flashes to have at least two sources. 

Once the valid meteorological flashes have been identified, we create gridded time- 
dependent 2-D horizontal maps of vertically-integrated flash origin density and flash extent 
density at 5 min intervals on a mesh that matches that of our WRF analyses. Flash origin 
density, which assigns unit value to the grid cell where a flash initiates, accurately tallies 
the actual flash rate in a storm, and is used here as a basis for calibration of WRF model 
proxies against forecast flash rates. Flash extent density, which assigns unit value to each 
grid cell a given flash traverses (Murphy and Demetriades 2005; Lojou and Cummins 2005), 
is better suited for gauging the spatial extent of the lightning threat, and is presented here 
as a means of comparing the areal coverage of observed and forecast lightning threat. It 
is also a field that is much less sparse than flash origin density, and thus more amenable 
to graphical display. For these reasons, it is used extensively in our figures describing our 
actual lightning observations. In all the LMA-based assessments that follow, we propose 
that a flash extent density of at least 1.0 flash traversal (5 min) -1 per grid box column 
represents an objective minimum threshold above which an observed lightning threat may 
be said to exist. 

Only the flash origin and extent density maps produced at those 5-min time intervals 
that coincide temporally with the WRF output analysis times are used in direct compar- 
isons with WRF data. These 5-min snapshots of lightning data are similar to what have 
been used by other researchers in their observational studies of lightning in storms, and 
represent a good compromise between the brevity needed to reveal trends in storm flash 
rate, and a duration long enough to provide adequate sampling of the lightning flash rate 
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in storms, especially in low flash-rate situations. 


c. Calibration of WRF lightning proxies 

Based on previous published work (e.g., Petersen et ah 2005), we propose that one 
useful estimate of the total flash rate may be based on the resolved upward flux wq g of 
large precipitating ice (i.e., graupel) in the mixed phase region at -15°C. We designate this 
first type of threat estimate by the symbol F\. For this threat we thus assume: 

-Fl = f{{m g )m) (i) 

where w is the vertical velocity, q g is the graupel mixing ratio, and the subscript m attached 
to the flux implies evaluation at the -15° C level in the mixed phase region. Based on 
published findings and our own results described below in section 3, we will show that 
the functional dependence symbolized by / in (1) is well described by a simple linear 
relationship. 

Because the WRF model, like many other cloud-resolving models, cannot be expected 
to simulate individual convective cells in exactly the same locations and times as observed 
storms, one cannot use simple scatterplots of individual simulated and observed grid point 
data to assess the nature of the functional dependence /. Instead, we assume that in 
successful simulations the model produces convection which, at its most intense, has a 
character that is a reasonably faithful rendition of the most intense storms observed any- 
where in the domain at any time during the simulation period. We thus compare the 
maximum gridded observed flash rate to the maximum simulated graupel flux within the 
simulation time frame to assess the functional form of /. The results of this comparison 
and calibration effort are discussed in more detail in Section 3a. 

The second model proxy field considered is based on the gridded vertical integral of 


13 



graupel, snow and cloud ice. Although the cloud ice field makes relatively little contribution 
to the amplitude of the lightning flash rate threat, it serves to distribute the associated 
lightning threat over a wider area, and thus may help account for extended anvil lightning 
events not readily captured by the graupel flux-based methods, which tend to be more 
confined to active updraft regions. Thus for our second threat F 2 we have: 

F 2 = kJ P(lg + Qs + qi)dz) (2) 

where p is local air density, and q g , q s and are the simulated mixing ratios of graupel, 
snow and ice, respectively, and h is a functional relationship that will be shown in Section 
3a to be a simple linear proportionality. 

3. Results 

Here we present results of our calibration efforts, followed by data from the two light- 
ning threat algorithms applied to two successful WRF simulation cases. One case is from 
00-06 UTC 30 March 2002, when an isolated tornadic supercell erupted in north Alabama 
equatorward of a broken cold frontal squall line moving through Tennessee. The supercell 
storm produced the highest peak lightning flash rate observed in any storm from this case, 
roughly 60 flashes min -1 just before generating an F-l tornado near Albertville, Alabama. 
Expressed in terms of flash origin density, the supercell produced a peak value of 12 flashes 
(5 min) -1 per grid box column. 

The second case is from 12-20 UTC 10 December 2004, when several small storms 
produced 2.5 cm hail after an intrusion of destabilizing, cold midtropospheric air advected 
across the Tennessee Valley region following a cold front. The flash rates from the strongest 
of these hailstorms were less than 5 flashes min -1 , and flash densities never exceeded 3 
flashes (5 min) -1 per grid box column. These two cases are chosen to span the approximate 
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range of severe convective storm event types found in the North Alabama region. Full 
LMA data were also available for both cases, as were the data needed to launch the WRF 
simulations. 

a. Calibration of threat indices 

Before interpreting results from our two selected case studies, it is necessary to com- 
pare findings from our successful simulations in order to establish reasonable calibration 
curves for the two lightning threat indices described herein. For threat Fi, which is based on 
WRF-simulated graupel fluxes at -15° C, we searched each simulation and its corresponding 
LMA-derived lightning flash rate observations, and extracted the maximum WRF graupel 
flux and LMA flash rate density from each, plotting them as single points in Fig. 1. Objec- 
tive least-squares fits of the data using reduced major axis regression techniques indicate 
that, for both Fi and F 2 , the 95% confidence limits contain the origin. Based on this, 
we devised linear regressions with intercepts at the origin for both our threats, and both 
regressions are found to be significant at more than the 99% confidence level. The slope 
thereby obtained for the Fi line is found to be 0.042. The raw cross-correlation between 
LMA flash rate and graupel flux is 0.67, with rather tight clustering of the data near lower 
values of each variable, but with enhanced scatter at large values. 

Based on the findings from Fig. 1, we propose that the functional relationship in (1) 
can be well described by the simple linear equation: 

F 1 = k 1 (wq g ) rn (3) 

where the coefficient k\ is the 0.042 slope of the line plotted on the figure. 

Following the calibration of lightning threat amplitude as in (3), we now consider the 
separate issue of threat areal coverage. To generate a threat field that has both reasonable 
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areal coverage and peak flash rate threat values, it is necessary to set a nonzero lower limit 
for the proxy threat field, which is equivalent to establishing a threshold contour for the 
field when plotted on a map. Values of the field that lie below this threshold are considered 
to pose no lightning threat, and are not contoured on our field map plots presented below. 
For both lightning threats, the optimum threshold contour levels are different, and are 
identified empirically via experimentation. For threat we find a threshold value of 
0.01 predicted flash origin points (5 min) -1 per grid box column to be the best that can 
reasonably be obtained. This threshold often leads to an underprediction of threat area, 
an issue which is discussed in more detail in section 5. However, it is probably as small a 
value as should be considered; sensitivity tests show that further reduction provides little 
or no increase in areal coverage of this threat. 

We follow a similar procedure in Fig. 2 for threat F 2 , which is based on gridded 
vertically integrated ice content from the WRF. Calibration procedures similar to those 
described above again indicate that the function h can be satisfactorily modeled as a 
linear proportionality, similar to what was found for threat F\ except for the value of the 
proportionality constant. The data points used in this calibration are shown in Fig. 2. 
For threat F 2 , we thus have: 

F 2 = k 2 J p(q g + q s + qi)dz (4) 

where the integration is over the full storm depth. From Fig. 2, our zero-intercept re- 
gression calculations indicate that the constant k 2 may be estimated as 0.20. For Fig. 
2, the objective raw cross-correlation between the LMA flash density and the values of 
vertically integrated ice is 0.83. In evaluating the areal coverage of threat F 2 , we find that 
a contouring threshold of 0.40 predicted flash origin points (5 min) -1 per grid box column 
provides approximately correct areal coverage for this threat. This threshold implies a 
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need for at least 2.1 kg m -2 of vertically integrated ice through the depth of storm anvils 
for the existence of a lightning threat there. 

b. Spring supercell and squall line 

At 0000 UTC on 30 March 2002, a warm, moist and unstable airmass was in place 
over the Tennessee Valley region, with a cold front approaching from the northwest. A 
capping inversion prevented storm formation until near sunset, when an isolated supercell 
erupted near Tupelo, Mississippi, and moved east into Alabama. Shortly afterwards, new 
severe storms quickly developed near the Smoky Mountains in eastern Tennessee, while 
other severe storms built along the cold front in middle Tennessee. The supercell passed 
south of Huntsville, Alabama, and went on to produce an F-l tornado west of Albertville, 
Alabama, after 0520 UTC. 

The WRF model produces a representative depiction of these developments during 
the 6-h forecast launched at 0000 UTC 30 March. A plot of the WRF sounding valid 
at (latitude, longitude) = (34.4, -88.1) and 0300 UTC 30 March 2002 (Fig. 3) indicates 
that CAPE of more than 2000 J kg -1 , along with moderate veering shear, was present 
in northern Alabama to support the severe convection. The reflectivity field (Fig. 4, 
gray shades) observed by the NWS Doppler radar at Hytop, Alabama, (KHTX) at 0400 
UTC indicates peak low-level reflectivities near or above 60 dBZ in both the supercell at 
(latitude, longitude) = (34.5, -86.7) and in portions of the other convective lines. LMA- 
derived flash origin densities reach a peak of 10 flashes (5 min) -1 per grid box column in 
the supercell at 0400 UTC, with an absolute peak of 12 having been reached earlier at 
0145 UTC. However, as stated earlier, because of the sparseness of the field of flash origin 
density, we opt to depict the lightning activity in terms of flash extent density (Fig. 4, 
contours). We also show on Fig. 4 the areal coverage fraction of the LMA-derived flash 
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extent density field exceeding our minimum plotted contour, which represents 1.0 flash 
traversal (5 mi n)" 1 per grid box column. The coverage fraction for this contour threshold 
turns out to be 0.148, which means that LMA observations indicate that 14.8% of the 
area of the analysis domain experienced actual lightning threat, either at the ground or 
overhead, over the 5-min period beginning at 0400 UTC. 

In Fig. 5 we present the fields of WRF-derived reflectivity at -15°C (grayshades) and 
predicted flash origin density F\ (color contours), valid at 0400 UTC 30 March. The areal 
coverage of the predicted F\ flash rate field in Fig. 6 is 0.085, which is approximately 60% 
of the actual observed flash extent density coverage in Fig. 4. It is, however, much larger 
than the coverage of the observed flash origin density for this event (not shown) , which is 
less than 0.005. On the other hand, the domainwide areal fractions of the predicted and 
observed lightning threats are much smaller than the areal fraction of positive CAPE, which 
is 0.64. In section 4, we will discuss the significance of the areal coverage discrepancies 
between the predicted and observed flash density. We emphasize here that our calibration 
methods yield predicted flash densities that are strictly referenced to observations of flash 
origin density, with the threat areal coverage referenced to the areal coverage of observed 
flash extent density. 

Fig. 6 contains a map of contoured lightning threat F 2 , based on the vertical integral 
of simulated ice, superimposed on the simulated field of anvil-level (near 9 km) cloud ice 
(shaded). Areal coverage fraction of this threat field, just under 0.158, is also indicated on 
the plot. It is somewhat larger than that for the first threat field, and is very close to the 
observed threat coverage shown in Fig. 4. 

To check for overall agreement of the intensity forecasts of the lightning threat, we 
present in Fig. 7 the time series of domain- wide peak values of LMA-derived flash origin 
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density along with the corresponding predicted peak values from the two threat methods 
Fi and FY From Fig. 7 we observe that WRF takes at least 60 min to begin to produce 
its first predicted lightning, with a large increase in predicted lightning at 90 min, which 
lags the observations in this regard. Once deep convection is fully established in the WRF 
model, both threat methods deliver flash rate density values within the range of the LMA 
observations, although the actual peak of LMA flash rate density occurs rather early on in 
the simulation period (just before 120 min), when WRF is still developing deep convective 
cells. There is substantial temporal variability in both the LMA-derived peak flash density 
and the predicted peak values of F i, but much less temporal variability in the predicted 
peak values of iY, which tend to plateau near their maximum values. 

Fig. 8 contains the time series plots of LMA-observed flash extent density areal 
coverage fraction that result when we apply the threshold of 1.0 flash traversals (5 min) -1 
per grid box column, our postulated lower limit in section 2b for observable lightning threat. 
Also plotted are the areal coverage fractions for threats F 1 and F 2 , each thresholded by 
the amounts mentioned in section 3a. All the areal coverage time series are more slowly 
varying than the peak threat value time series of Fig. 7. Of additional interest is the 
extended time lag seen in the onset of significant areal coverage of lightning threat, as 
compared to the onset times of the intense but localized lightning threat peaks in Fig. 7. 
There are also more general temporal mismatches in the times of local peak observed and 
predicted lightning threat coverage, reinforcing our earlier caveats about the limitations 
of single individual numerical simulations as forecast tools for real-world convective-scale 
events. Finally, the peak areal coverage values for threat F 1 are still somewhat too small 
compared to observed coverage values, suggesting that the threshold contours arrived at 
for F 1 (see Section 3a) might be too large. However, because this proxy threat field is so 
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closely tied to strong updraft cores, it was found that lowering its threshold failed to yield 
useful increases in its corresponding areal coverage. 

c. Cold season hailstorm 

On 10 December 2004, a cold airmass aloft moved into the Tennessee Valley, with 
sufficient lingering warmth and moisture at low levels to allow for the formation of shallow 
but moderately intense convection. The storms formed just after midmorning (1530 UTC) 
in Tennessee, with additional storms forming in north Alabama around local noon (1800 
UTC). With a cold atmospheric profile aloft (Fig. 9), the storms had approximately 500 
J kg -1 of CAPE with which to build updrafts. Some of the storms produced hail at the 
surface as large as 2.4 cm in diameter, according to reports in Storm Data. Maximum 
low-level reflectivities observed by KHTX reached approximately 58 dBZ (see Fig. 10 for 
representative data from 1900 UTC). The total lightning flash rates from the several storms 
remained rather small, however, never exceeding 5 flashes (min) -1 even in the strongest 
cells. Observed peak flash origin densities never exceeded 3 flashes (5 min) -1 per grid box 
column, as suggested by the modest amplitudes of the flash extent density field. Observed 
flash extent density in excess of 1.0 traversals (5 min) -1 per grid box displays areal coverage 
of only 0.018. Both the peak flash densities and areal coverages for this event are more 
than a factor of 4 smaller than those for the 30 March 2002 case. 

In the 8-h simulation launched at 1200 UTC 10 December 2004, the WRF model does 
well at capturing the location, intensity and character of the deep convection, but shows 
some error in the timing of its initiation. In the WRF results, for example, deep convection 
begins after t — 200 min (around 1520 UTC) in Tennessee. In the observations, however, 
a few brief, weak storms are seen there as early as t = 60 min (1300 UTC), becoming 
more widespread only after t = 320 min (1720 UTC). By t = 360 min (1800 UTC), both 
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the WRF and observations indicate additional deep convection developing in northeastern 
Alabama. All the storms occur in the form of small isolated cells, with some clustering into 
short lines. At 1900 UTC, the field of surface-based pseudoadiabatic CAPE from WRF 
(not shown) reveals a large region of values reaching near 400 J kg -1 , with a predicted 
lightning threat field F\ concentrated into small patches associated with just a few of 
the cells (Fig. 11). As with the 30 March 2002 case, areal coverage of the predicted 
lightning flash density F\ (0.012) falls somewhat short of the observed coverage of flash 
extent density (0.018). In addition, the predicted and observed lightning threat areas are 
again much smaller than the area of positive CAPE (domain fractional coverage of 0.60) 
simulated by WRF. 

Fig. 12 shows the anvil- level cloud ice field near z — 5 km, along with the computed 
lightning threat field F 2 , based on column-integrated ice. As in Fig. 6, areal coverage of the 
threat is expanded relative to the other two threat fields, reaching a fractional coverage of 
0.016. Even though there is considerable ice aloft in the western half of the domain, threat 
F 2 correctly predicts that all the lightning is confined to the northeastern and eastern 
portions of the WRF simulation region. 

The 8-h time series of domainwide maximum observed flash origin density for the 
10 December 2004 event is presented in Fig. 13, along with the time series of maximum 
predicted flash origin density from methods Fi (THREAT1) and F 2 (THREAT2). The 
time series plots reiterate that WRF is approximately 2 h too fast in developing multiple 
deep convective storms in the domain, although a few isolated weak observed storms are 
found even before the onset of the multiple simulated storms. The peak flash rate densities 
of threat F\ are only about half those of the observations, while threat F 2 yields peak 
flash rate densities much closer to those observed. None of the predicted peak flash rates 


21 



from Figs. 7 and 13, however, matches the observed peaks exactly, because of the fact 
that neither of the two cases presented here have predicted peak values that fall directly 
on the consensus calibration curves of Figs. 1-2. 

Fig. 14 portrays the 8-h time series of the areal coverage fractions of observed and 
predicted lightning threat for this case. Note that the vertical axis of Fig. 14 has a 
much smaller range than Fig. 8, because of the much smaller and more localized lightning 
threat coverage in the 10 December 2004 event. The figure shows clearly that the weak, 
early storms seen in the observations have negligible areal coverage in the domain, and 
that the WRF model develops more widespread storminess approximately 2 h earlier than 
indicated by the observations. The limited areal coverage of lightning threat in this event 
also poses a challenge to our prognoses of threat coverage, with Threat F 2 exhibiting a 
fractional peak of 0.023, some 26% smaller than the observed peak of 0.031. As with the 
much more vigorous 30 March 2002 case, Threat F\ predicts coverage even smaller than 
F 2 . As before, however, both our prognostic threat fields exhibit much more restricted 
and specific coverage than more traditional convection indicators such as positive CAPE, 
whose coverage fraction at 19 UTC, the time of Figs. 10-12, is 0.60. 

4. Blending of calibrated threats 

Results from the previous section demonstrate that threat F±, which is based on 
WRF graupel flux at -15°C, captures much of the temporal variability of observed peak 
lightning flash density in the mature convection, but underestimates the areal coverage 
of the actual lightning activity. Threat F 2 , meanwhile, based on vertically integrated ice, 
including anvil cirrus, captures the areal coverage of lightning activity well, but portrays 
the temporal variability of the lightning flash rates less accurately than does F 1 , because 
it is a vertically integrated quantity. We are thus motivated to construct a blended threat 
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index that retains the temporal sensitivity of F\ and better captures the areal coverage of 

f 2 . 

Because both F\ and jF 2 have been calibrated against peak lightning flash rate density, 
any weighted combination of the two threat fields will also be properly calibrated as well. 
To retain as much of the temporal variability inherent in threat F \ , however, any blended 
threat should probably be weighted heavily toward F \ . On the other hand, even a lightly 
weighted contribution by the F 2 field should suffice to provide the desired increase in net 
areal coverage to the blended threat F 3 , after proper thresholding. We thus propose that 
a workable blended threat F% could be based upon: 

F 3 = r\Fi + r 2 F 2 (5) 

where r\ = 0.95 and r 2 = 0.05, based on results of sensitivity tests of the effects of various 
weight choices on resulting peak flash rate densities and areal coverages. Application of 
(5) to the WRF data for 04 UTC, 30 March 2002, yields the threat field F 3 shown in Fig. 
15. Application of a minimum threshold flash rate density of 0.02 flash origins (5 min) -1 
per grid box column - the first contour plotted in Fig. 15 - yields very good agreement 
between the predicted (0.17) and observed (0.15) lightning threat coverages. The time 
series of predicted and observed peak flash densities for F 3 for the 30 March 2002 case are 
shown in Fig. 16, while the predicted and observed threat areal coverage time series are 
given in Fig. 17. As expected and desired, the areal coverage for jF 3 is larger than for 
Fi (see Fig. 5), and approaches that of F 2 , while the accuracy of the F 3 peak flash rate 
density is unchanged, even as most of the temporal sensitivity of Fi is retained (see Fig. 
17). 

The blended threat F 3 works almost as well for the low flash-rate case of 10 December 
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2004. Fig. 18 shows the field of predicted threat F3 at 19 UTC, while Figs. 19 and 20 
present the time series of observed and predicted peak flash rate densities and threat areal 
coverages, respectively. The same threshold as in Figs. 15 and 17, 0.02 flash origins (5 
min) -1 per grid box column, was used in assessing threat areal coverage. Although F3 
still underpredicts the peak areal coverage for this case (0.02 versus 0.03 observed), the 
threat predictions based on F3 perform better than those from F\ or F2 alone, despite the 
inherent challenges posed by the very low flash rates. 

5. Discussion 

In general, caution is always warranted when dealing with small sample sizes, as is 
the case in the calibration efforts herein. In both our calibration charts (Figs. 1 and 2), we 
have modeled the relationship between the two variables as linear, but it is possible that 
nonlinear behavior might emerge if a larger data sample were available. This is considered 
a possibility in light of the behavior of the rightmost data point in Fig. 1, which has 
a graupel flux of 400 m s -1 , but falls below the regression line. This particular point 
derives from a simulation for a severe squall line that occurred on 31 May 2004, and is 
characterized by our largest updraft strengths and largest environmental CAPE values. It 
is possible that the model’s lack of representation of the hail species may have led to an 
overproduction of graupel for this intense storm system simulation, although additional 
cases featuring very large CAPE would be needed to confirm this speculation. 

Despite the uncertainties attending our small sample of cases, the methods used to 
devise the lightning threat proxy fields F 1 and F 2 still appear to yield useful approximate 
indicators of the lightning threat. Our blended threat F 3 appears to mitigate the short- 
comings of the threats F\ and F 2 , by using a judicious weighted average of the two basic 
threat fields. By design, F 3 tends to match the peak flash rate densities in storm systems, 
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but it also provides a realistic match for lightning threat areal coverage, while retaining 
most of the temporal variability of lightning in observed mature storms. It should be noted, 
however, that our threat F 3 should not be expected to match either the observed peak 
flash rate densities or threat areal coverages exactly for any specific storm cases, because 
of the statistical character of our calibration methods. For example, the peak values of F 3 
in Figs. 16 and 19 both lie below those of the observed peak flash densities. This is simply 
because both these observed cases happened to produce peak flash rate densities that fell 
above the least squares best fit calibration line of Fig. 1. The behavior of the threat 
areal coverage can differ, however, from that of peak flash rate density, as seen in Figs. 
17 and 20, where the predicted peak areal coverages are above and below those observed, 
respectively. When examining the values in instantaneous snapshots of data, matches can 
exhibit even more error, owing to the inability of the model to place convective storms of 
exactly the right intensity in exactly the right place at exactly the right time. 

All our threats also exhibit a much more confined areal coverage than other traditional 
environmental measures used in predicting the likelihood of thunderstorms, such as CAPE. 
Although the areal coverage of CAPE often overestimates the areal coverage of lightning 
threat, as in cases where capping or large-scale subsidence are present, it is often considered 
an important field to examine, because of the possibility that small-scale forcing might 
overcome the capping and allow for deep convection. For the cases presented here, positive 
CAPE existed over approximately 60% of the domain at any given time, but the actual 
instantaneous lightning threat area, as measured by the LMA flash extent density fields and 
reproduced by threat F 3, never exceeded 16% coverage. When viewed in a time- integral 
sense over the course of all seven simulations, CAPE in excess of 100 J kg -1 covered 88- 
100% of the domain, but time-integrated lightning threat covered only about 10-30% of the 
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domain for cases not involving major squall lines, and about 50-80% for squall line cases. 
These statistics give some idea of how important it is to avoid unnecessary overestimation 
of the area subject to the forecast lightning threat. 

Interestingly, our blended lightning threat coverage tended to be similar to the cover- 
age of the 20 dBZ and greater reflectivity area, and noticeably larger than the area within 
the 35 dBZ reflectivity isopleth. The latter, for example, was only 11% for the storms 
in Fig. 5, while the actual observed lightning threat coverage at that time was approx- 
imately 15%. The 35 dBZ area thus appears to resemble more closely the area covered 
by our graupel-flux-based threat F ± alone (8.5%), although the former is still somewhat 
larger than the latter. Nevertheless, our observations of lightning threat indicate that ice 
detrained in storm anvils, which is accounted for by threat F 2 and its contribution to the 
blended threat F 3 , should not be neglected in the assessment of total lightning threat. 

Note that although the blended lightning threat area is close to being coincident with 
the area of greater than 20 dBZ echo, the use of simple radar reflectivity thresholds to 
estimate lightning threat is far from straightforward. Our tests of reflectivity as a possible 
lightning proxy field suggest that nonlinearities are present in the calibration curves; in 
addition, there are difficulties in finding storm cases where the peak reflectivity did not 
migrate to high values, frustrating our attempts at constructing evenly populated calibra- 
tion curves. Likewise, it is difficult to compare our apparent dBZ-based lightning threat 
spatial extent thresholds with other research dealing with dBZ thresholds for lightning 
onset (see, e.g., Buechler and Goodman 1990), which commonly infer a need for at least 40 
dBZ in the mixed phase region. However, our finding that lightning threat areal coverage 
approximately matches the area inside the 20 dBZ reflectivity isopleth is not inconsistent 
with the 40 dBZ lightning onset threshold; an interpretation consistent with our data is 
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that, while threat areal coverage roughly matches the 20 dBZ or greater area, at least 

i 

some part of the threat zone must contain 40 dBZ or greater in the mixed phase region. 
Another way of viewing this situation is that lightning generally requires at least 40 dBZ 
for initiation, but it can often propagate out into regions of as little as 20 dBZ reflectivity. 

Our lightning threat computation methods also exhibit a positive bias in the mean, 
which is visible in the various time series plots of predicted and observed peak flash density 
(Figs. 7, 13, 16 and 19). This bias is an inevitable consequence of our choice to design 
our threats to match the peak flash densities, rather than the averages, in our storm cases. 
Although our main justification for this design choice is that the most intense storms 
are of greater concern than others, a secondary justification is that it is problematic to 
define a meaningful “average” flash density. Our methodology thus tends to overpredict 
mean flash densities, but this also provides a conservative design, from a warning point 
of view. Similarly, our threat methods will, when integrated spatially over a storm cell’s 
instantaneous footprint, tend to produce small overestimates of storm cell flash rate. 

The calibration of our prognosed fields of lightning threat is simple, but not without 
issues requiring careful attention. In particular, the determination of the calibration con- 
stants should be performed anew whenever a new model or significantly different version 
of a previously used model is employed. Previous work (Weisman et al. 1997; Bryan et 
al 2003) has shown how simulated field quantities can vary in amplitude as the model 
mesh is made coarser or finer. At this point, there is no general method for predicting 
quantitatively the amount of variation in field amplitudes that is realized when the model 
mesh changes. All that can be said is that, when comparing to quantities derived from 
the 2 km native model mesh used here, most field quantities will grow (shrink) in peak 
amplitude, probably by some tens of percent, as the model mesh is, say, halved (doubled). 
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We thus emphasize that the specific calibration constants presented herein may not gen- 
erally provide optimum results for model configurations other than the one discussed in 
this study. On the other hand, the general approaches that we have outlined here are 
firmly based on storm physics and should have some general regime-independent validity. 
As cloud models, initialization procedures, and data assimilation techniques improve, and 
more lightning cases are archived and analyzed, we expect that the functional relations 
and diagnosed constants described here will evolve toward greater accuracy and generality, 
and that overall model capability to predict lightning threat quantitatively will improve. 

Although we found that the WRF model did not always produce successful simulations 
of observed convection, even for some severe convective events, the difficulties in obtaining 
reliable simulations of less intense convective events were even more challenging. Yet both 
high flash-rate and low flash-rate cases are needed in order to make robust algorithm 
calibrations. We continue to explore our data archives for other low-end cases that might 
yield additional insights. L 

6. Summary and Outlook 

We have presented evidence that application of simple physics-based lightning con- 
cepts to the calibration of fields output by cloud-resolving models can provide quantita- 
tively calibrated maps of lightning threat for use by forecasters. These simple methods 
have the potential to allow for useful short-term lightning threat forecasts, without the 
need for adding expensive and complex electrification subroutines to cloud-resolving mod- 
els. Drawbacks include the fact that individual simulations tend to do an imperfect job 
at locating and timing convective storms, and that the cloud-resolving model output that 
is used to forecast lightning threat may need to be recalibrated against observed lightning 
data for major changes to the model grid mesh and physics. Intrinsic model errors also 
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introduce uncertainty into the values of the calibration constants needed in constructing 
model-based lightning threat forecasts. 

Many of the shortcomings of the present WRF model are likely to be alleviated, at least 
in part, by future improvements to the model physical parameterizations. In particular, a 
more sophisticated double-moment cloud microphysics scheme that includes at least some 
form of hail in addition to graupel, would be desirable. Of course, if more species of large 
precipitating ice are included, our algorithms must be modified to include their effects, 
even though the graupel species, along with small hail, likely do most of the work in the 
charge separation process. Running the model on a finer grid mesh would probably also be 
beneficial in removing the small weak bias seen in the current model storms. Even when 
these improvements are made, however, some significant errors may well remain, owing 
to the imperfections in the initial and lateral boundary condition fields used to run the 
simulations. 

In the future, we suggest that some of the uncertainties associated with the model 
data can be addressed by running small ensembles using diverse choices for microphysics 
and boundary layer schemes, and, if possible, initial conditions. Probability-based analysis 
of the patterns of convection found across such ensembles will allow for objective statistical 
assessment of the skill of our lightning threat methods in prospective future case studies 
or operational settings. Although studies using results from ensembles of simulations 
are considered a top priority, experiments with newer versions of the WRF model and 
other cloud-resolving models are also contemplated. In addition, we seek to identify and 
analyze other new cases of very low, intermediate, and very high flash-rate storm events, 
and to run cloud-resolving simulations on them for the purpose of refining our empirical 
estimates of the functional forms and curve-fitting coefficients of our various lightning 
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threat predictors. In particular, we are also in the process of reanalyzing the database 
of Cecil et al. (2005) to obtain lightning-reflectivity relationships based on gridded data, 
rather than storm system precipitation features, and to construct additional lightning 
threats based on reflectivity profile shape. Finally, the existence of simple, apparently 
invertible linear relationships between total flash rate density and parameters such as 
graupel flux and vertically integrated ice suggest opportunities for ways to convert observed 
lightning flash rate data from satellite-based sensors such as GLM into forms that can be 
assimilated into operational forecast models. 
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TABLE 1 


Model numerical and physical parameters 


Parameter 

Horizontal resolution Ax, Ay 
Number of vertical sigma levels 
Large time step At 
Dynamical core , 

PBL turbulence scheme 
Shortwave radiation scheme 
Longwave radiation scheme 
Land surface model scheme 
Microphysics scheme 
Graupel density 
Basic initialization fields 
Extra initialization data used 
Lateral boundary conditions 


Value 
2000 m 
51 

12.0 s 

Eulerian mass 
YSU 
Dudhia 
RRTM 
Noah 

WSM6, graupel 
300 kg m~ 3 

AWIP212 NCEP EDAS 
METAR,ACARS,WSR88D 
Eta 3-h forecasts 


YSU = Yonsei University Scheme 
RRTM = Rapid Radiative Transfer Model 

METAR = regular hourly Surface aviation weather observations 


ACARS = Aircraft Communications Addressing and Reporting System 



FIGURE CAPTIONS 


FIG. 1. Calibration scatterplot comparing domainwide maxima for observed LMA flash 
rate density (y-axis; units are flashes (5 min) -1 per grid box column) and WRF-simulated 
proxy data for lightning threat method 1 (x-axis; units are ms” 1 ), based on WRF-derived 
graupel flux at the -15° C level. A straight line passing through the origin and having a 
slope of 0.042 fits the data well. 

FIG. 2. Calibration scatterplot comparing domainwide maxima for observed LMA flash 
rate density (y-axis; units are flashes (5 min) -1 per grid box column) and WRF-simulated 
proxy data for lightning threat method 2 (x-axis, units are kg m -2 ), based on WRF-derived 
vertically integrated ice. A straight line passing through the origin and having a slope of 
0.20 fits the data well. 

FIG. 3. WRF skew-T log-p sounding at (latitude, longitude) = (34.4, -88.1), at 0300 UTC 
30 March 2002, showing surface-based pseudoadiabatic CAPE greater than 2000 J kg -1 
in the warm sector southwest of an observed tornadic supercell. Wind vectors are shown 
on the staff at right, with each full barb equal to 5 ms” 1 of windspeed. 

FIG. 4. Low-level radar reflectivity from KHTX Doppler radar at 0400 UTC 30 March 
2002 (grayshades) , and LMA-derived integrated flash extent density (contours) for a 5- 
min period at the same time. Flash extent density is depicted for clarity, although all 
calibrations and predictions of lightning threat are based on observations of flash origin 
density. 

FIG. 5. WRF-derived reflectivity at the -15° C level at 0400 UTC 30 March 2002 
(grayshades), and WRF-predicted flash origin density (contours) for a 5-min period at 
the same time, based on WRF graupel flux at the -15°C level. Instantaneous areal cover- 



age of predicted flash density is printed at the bottom of the figure. 

FIG. 6. WRF-derived anvil-level cloud ice field at 0400 UTC 30 March 2002 (grayshades) , 
and WRF-predicted flash origin density (contours) for a 5-min period at the same time, 
based on WRF vertically integrated ice mixing ratios. Instantaneous areal coverage of 
predicted flash density is printed at the bottom of the figure. 

FIG. 7. Time series plots of domainwide peak observed LMA flash origin density and 
predicted flash origin density based on two proxies: graupel flux (THREAT 1) and vertically 
integrated ice content (THREAT2), for the 30 March 2002 case. Time series data symbols 
are explained on the figure. 

FIG. 8. Time series plots of domainwide areal coverage fractions of LMA flash extent 
density and predicted flash origin density based on two proxies: graupel flux (THREAT1) 
and vertically integrated ice content (THREAT2), for the 30 March 2002 case. Time series 
data symbols are explained on the figure. 

FIG. 9. WRF skew-T log-p sounding at (latitude, longitude) = (34.8, -85.9), at 1800 UTC 
10 December 2004, showing surface-based pseudoadiabatic CAPE near 500 J kg -1 in the 
vicinity of an observed postfrontal hailstorm. Wind vectors are shown on the staff at right, 
with each full barb equal to 5 m s -1 of windspeed. 

FIG. 10. Low-level radar reflectivity from KHTX Doppler radar at 1900 UTC 10 December 
2004 (grayshades), and LMA-derived integrated flash extent density (contours) for a 5- 
min period at the same time. Flash extent density is depicted for clarity, although all 
calibrations and predictions of lightning threat are based on observations of flash origin 
density. 

FIG. 11. As in Fig. 5, but for 1900 UTC 10 December 2004 case. 



FIG. 12. As in Fig. 6, but for 1900 UTC 10 December 2004 case. 

FIG. 13. As in Fig. 7, but for 1200-2000 UTC 10 December 2004 case. 

FIG. 14. As in Fig. 8, but for 1200-2000 UTC 10 December 2004 case. 

FIG. 15. As in Fig. 5, but for blended threat (THREAT3), 0400 UTC 30 March 2002 
case. 

FIG. 16. As in Fig. 7, but for blended threat (THREAT3), 0000-0600 UTC 30 March 
2002 case. 

FIG. 17. As in Fig. 8, but for blended threat (THREAT3), 0000-0600 UTC 30 March 
2002 case. 

FIG. 18. As in Fig. 5, but for blended threat (THREAT3), 1900 UTC 10 December 2004 
case. 

FIG. 19. As in Fig. 7, but for blended threat (THREAT3), 1200-2000 UTC 10 December 
2004 case. 

FIG. 20. As in Fig. 8, but for blended threat (THREAT3), 1200-2000 UTC 10 December 


2004 case. 



Threat 1 Max, Graupel Flux at T = -15C 



FIG. 1. Calibration scatterplot comparing domainwide maxima for observed LMA flash 
rate density (y-axis; units are flashes (5 min) -1 per grid box column) and WRF-simulated 
proxy data for lightning threat method 1 (x-axis; units are ms -1 ), based on WRF-derived 
graupel flux at the -15°C level. A straight line passing through the origin and having a 
slope of 0.042 fits the data well. 



Threat2 Max, Vert Integrated Ice 



MAX VERT INTEGRATED ICE 


FIG. 2. Calibration scatterplot comparing domainwide maxima for observed LMA flash 
rate density (y-axis; units are flashes (5 min) -1 per grid box column) and WRF-simulated 
proxy data for lightning threat methods (x-axis, units are kg m -2 ), based on WRF-derived 
vertically integrated ice. A straight line passing through the origin and having a slope of 
0.20 fits the data well. 




FIG. 3. WRF skew-T log-p sounding at (latitude, longitude) = (34.4, -88.1), at 0300 UTC 
30 March 2002, showing surface-based pseudoadiabatic CAPE greater than 2000 J kg -1 
in the warm sector southwest of an observed tornadic supercell. Wind vectors are shown 
on the staff at right, with each full barb equal to 5 m s _1 of windspeed. 
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FIG. 4. Low-level radar reflectivity from KHTX Doppler radar at 0400 UTC 30 March 
2002 (grayshades) , and LMA-derived integrated flash extent density (contours) for a 5- 
min period at the same time. Flash extent density is depicted for clarity, although all 
calibrations and predictions of lightning threat are based on observations of flash origin 





FIG. 5. WRF-derived reflectivity at the -15°C level at 0400 UTC 30 March 2002 
(grayshades) , and WRF-predicted flash origin density (contours) for a 5-min period at 
the same time, based on WRF graupel flux at the -15° C level. Instantaneous areal cover- 
age of predicted flash density is printed at the bottom of the figure. 






FIG. 6. WRF-derived anvil-level cloud ice mixing ratio at 0400 UTC 30 March 2002 
(grayshades) , and WRF-predicted flash origin density (contours) for a 5-min period at the 
same time, based on WRF vertically integrated ice content. Instantaneous areal coverage 
of predicted flash density is printed at the bottom of the figure. 
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FIG. 7. Time series plots of domainwide peak observed LMA flash origin density and 
predicted flash origin density based on two proxies: graupel flux (THREAT 1) and vertically 
integrated ice content (THREAT2), for the 30 March 2002 case. Time series data symbols 
are explained on the figure. 


LMA, WRF FLASH AREA FRAC vs time, 20020330 



FIG. 8. Time series plots of domainwide areal coverage fractions of LMA flash extent 
density and predicted flash origin density based on two proxies: graupel flux (THREAT1) 
and vertically integrated ice content (THREAT2), for the 30 March 2002 case. Time series 
data symbols are explained on the figure. 
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FIG. 10. Low-level radar reflectivity from KHTX Doppler radar at 1900 UTC 10 December 
2004 (grayshades), and LMA-derived integrated flash extent density (contours) for a 5- 
min period at the same time. Flash extent density is depicted for clarity, although all 
calibrations and predictions of lightning threat are based on observations of flash origin 
density. 


FIG. 11. WRF-derived reflectivity at the -15°C level at 1900 UTC 10 December 2004 
(grayshades) , and WRF-predicted flash origin density (contours) for a 5-min period at the 
same time, based on WRF graupel flux at the -15°C level. Instantaneous areal coverage 
of predicted flash density is printed at the bottom of the figure. 
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FIG. 12. WRF-derived anvil-level cloud ice mixing ratio at 1900 UTC 10 December 2004 
(grayshades) , and WRF-predicted flash origin density (contours) for a 5-min period at the 
same time, based on WRF vertically integrated ice content. Instantaneous areal coverage 
of predicted flash density is printed at the bottom of the figure. 
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FIG. 13. Time series plots of domainwide peak observed LMA flash origin density and 
predicted flash origin density based on two proxies: graupel flux (THREAT 1) and vertically 
integrated ice content (THREAT2), for the 10 December 2004 case. Time series data 
symbols are explained on the figure. 
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FIG. 14. Time series plots of domainwide areal coverage fractions of LMA flash extent 
density and predicted flash origin density based on two proxies: graupel flux (THREAT1) 
and vertically integrated ice content (THREAT2), for the 10 December 2004 case. Time 
series data symbols are explained on the figure. 
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FIG. 16. As in Fig. 7, but for blended threat (THREAT3), 0000-0600 UTC 30 March 
2002 case. 
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FIG. 17. As in Fig. 8, but for blended threat (THREAT3), 0000-0600 UTC 30 March 
2002 case. 










FIG. 19. As in Fig. 7, but for blended threat (THREAT3), 1200-2000 UTC 10 December 
2004 case. 



FIG. 20. As in Fig. 8, but for blended threat (THREAT3), 1200-2000 UTC 10 December 
2004 case. 






